﻿#pragma execution_character_set("utf-8")
#include <stdexcept>
#include "TheTask.h"
#include "gdal_priv.h"
#include <vector>
using std::vector;
using std::runtime_error;
#include "modisutitls.h"
#include <qfile.h>
#include <fstream>
#include "wdirtools.h"

#define THETASK_NODATA -9999

TheTask::TheTask(string infilename, string outfilename,
	double a1, double a2, double a3, double a4, bool useWarp, 
	string gdalwarp, string vrtTemFilename)
{
	_infilename = infilename;
	_outfilename = outfilename;
	_a1 = a1;
	_a2 = a2;
	_a3 = a3;
	_a4 = a4;
	_useWarp = useWarp;
	_gdalwarp = gdalwarp;
	_vrtTemFilename = vrtTemFilename;
}
void TheTask::run()
{
	GDALDataset* emmds = 0;
	GDALDataset* szads = 0;
	GDALDataset* outds = 0;
	try {
		string emmdsname = string("HDF4_EOS:EOS_SWATH:") + _infilename + ":MODIS_SWATH_Type_L1B:EV_1KM_Emissive";
		string szadsname = string("HDF4_EOS:EOS_SWATH:") + _infilename + ":MODIS_SWATH_Type_L1B:SolarZenith";
		vector<double> scaleArr;
		vector<double> offsetArr;
		modisUtils_getRadianceScaleAndOffset(emmdsname, scaleArr, offsetArr);
		if (scaleArr.size() != 16 || offsetArr.size() != 16) {
			throw runtime_error("either size of scaleArr nor offsetArr is not 16.");
		}

		string error;
		emmds = (GDALDataset*)GDALOpen(emmdsname.c_str(), GA_ReadOnly);
		if (emmds == 0) {
			throw runtime_error("打开输入辐亮度数据集失败。");
		}
		szads = (GDALDataset*)GDALOpen(szadsname.c_str(), GA_ReadOnly);
		if (szads == 0) {
			throw runtime_error("打开输入天顶角数据集失败。");
		}
		int xsize = emmds->GetRasterXSize();
		int ysize = emmds->GetRasterYSize();
		int numbands = emmds->GetRasterCount();
		if (numbands != 16) {
			throw runtime_error("输入辐亮度数据集波段数不是16，不是标准MODIS产品文件。");
		}

		string outfilename1 = _outfilename;
		if (_useWarp) {
			outfilename1 += ".tmp.tif";
		}
		GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
		outds = driver->Create(outfilename1.c_str(), xsize, ysize, 1, GDT_Float32, 0);

		//sza
		int xsize1 = szads->GetRasterXSize();
		int ysize1 = szads->GetRasterYSize();
		vector<float> szadata(xsize1*ysize1);
		szads->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize1, ysize1, szadata.data(), xsize1, ysize1, GDT_Float32, 0, 0, 0);
		float x_1km_to_5km = 0.2f ;
		float y_1km_to_5km = 0.2f;

		//processing
		int per0 = -1;
		double perScale = 100.0 / ysize;
		vector<unsigned short> raddata31(xsize);
		vector<unsigned short> raddata32(xsize);
		
		vector<float> outdata(xsize);
		const double C1 = 1.191042E8;
		const double C2 = 1.4387752E4;
		const double lamad31 = 11.03;
		const double lamad32 = 12.02;
		const double lamad31_pow5 = lamad31 * lamad31 * lamad31 * lamad31 * lamad31;
		const double lamad32_pow5 = lamad32 * lamad32 * lamad32 * lamad32 * lamad32;
		const double DEG2RAD = 3.1415926 / 180.0;
		for (int iy = 0; iy < ysize; ++iy) {
			emmds->GetRasterBand(11)->RasterIO(GF_Read, 0, iy, xsize, 1, raddata31.data(), xsize, 1, GDT_UInt16, 0, 0, 0);
			emmds->GetRasterBand(12)->RasterIO(GF_Read, 0, iy, xsize, 1, raddata32.data(), xsize, 1, GDT_UInt16, 0, 0, 0);
			for (int ix = 0; ix < xsize; ++ix) {
				int si31 = raddata31[ix];
				int si32 = raddata32[ix];
				int szax = ix * x_1km_to_5km;
				int szay = iy * y_1km_to_5km;
				float szaRad = szadata[szay * xsize1 + szax]* (0.01 * DEG2RAD );
				float outvalue = THETASK_NODATA;
				if (si31 > 0 && si31 < 32700 && si32>0 && si32 < 32700) {
					double r31 = scaleArr[10] * (si31 - offsetArr[10]);
					double r32 = scaleArr[11] * (si32 - offsetArr[11]);
					//计算亮温 https://www.sciencedirect.com/science/article/pii/S0303243418306287 波长波数取值来源 wv31=11.03um wv32=12.02um
					// C1 C2 from https://ncc.nesdis.noaa.gov/Planck.php
					double temp31 = C2 / (lamad31 * log(C1 / lamad31_pow5 / r31 + 1.0)); // unit K
					double temp32 = C2 / (lamad32 * log(C1 / lamad32_pow5 / r32 + 1.0));
					//WST = A1 + A2*T31 + A3*|T31-T32| + A4*[SEC(t0)-1]*(T31-T32)
					double wst = _a1 + _a2 * temp31 + _a3 * fabs(temp31 - temp32) + _a4 * (1.0 / cos(szaRad) - 1.0) * (temp31 - temp32);
					outvalue = wst;
				}
				outdata[ix] = outvalue;
			}
			outds->GetRasterBand(1)->RasterIO(GF_Write, 0, iy, xsize, 1, outdata.data(), xsize, 1, GDT_Float32, 0, 0, 0);
			int per1 = iy * perScale;
			if (per1 != per0) {
				per0 = per1;
				this->notifyProgress(per0, 100);
			}
		}

		outds->GetRasterBand(1)->SetNoDataValue(THETASK_NODATA);
		if (emmds) GDALClose(emmds); emmds = 0;
		if (szads) GDALClose(szads); szads = 0;
		if (outds) GDALClose(outds); outds = 0;

		if (_useWarp) {
			string warpError;
			if (warp(outfilename1,xsize,ysize,warpError)) {
				this->notifyProgress(100, 100);
				this->notifyDone();
			}
			else {
				string error1 = string("坐标变换失败：") + warpError;
				this->notifyError(error1);
			}
		}
		else {
			this->notifyProgress(100, 100);
			this->notifyDone();
		}
	}
	catch (runtime_error& ex) {
		if (emmds) GDALClose(emmds); emmds = 0;
		if (szads) GDALClose(szads); szads = 0;
		if (outds) GDALClose(outds); outds = 0;
		string error = ex.what();
		this->notifyError(error);
	}
}

bool TheTask::warp(const string tifFilename, const int wid, const int hei, string& error)
{
	string temtext;
	bool tok = wDirTools::readFullText(_vrtTemFilename, temtext);
	if (tok == false) {
		error = "读取VRT模板文件失败。";
		return false;
	}
	QString vrtContent = QString::fromStdString(temtext) ;
	vrtContent = vrtContent.replace("$WID", QString::number(wid));
	vrtContent = vrtContent.replace("$HEI", QString::number(hei));
	vrtContent = vrtContent.replace("$HDF", _infilename.c_str());
	vrtContent = vrtContent.replace("$TIF", tifFilename.c_str());
	string content = vrtContent.toStdString();

	string vrtfilename =  tifFilename  + ".vrt";
	std::ofstream ofs(vrtfilename);
	if (ofs.good() == false) {
		error = "创建VRT文件失败。";
		return false;
	}
	ofs << vrtContent.toStdString();
	ofs.close();

	string cmd1 = _gdalwarp + " -overwrite -tr 0.01 0.01 -t_srs WGS84 " + vrtfilename + " " + _outfilename;
	system(cmd1.c_str());

	wDirTools::removeFile(vrtfilename);
	wDirTools::removeFile(tifFilename);
	emit onTaskDone();

	return true;
}